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ABSTRACT 

We study the effects of anisotropic thermal conduction along magnetic field lines on 
an accelerated contact discontinuity in a weakly collisional plasma. We first perform a 
linear stability analysis similar to that used to derive the Rayleigh- Taylor instability 
(RTI) dispersion relation. We find that anisotropic conduction is only important for 
compressible modes, as incompressible modes are isothermal. Modes grow faster in the 
presence of anisotropic conduction, but growth rates do not change by more than a 
factor of order unity. We next run fully non-linear numerical simulations of a contact 
discontinuity with anisotropic conduction. The non-linear evolution can be thought of 
as a superposition of three physical effects: temperature diffusion due to vertical con- 
duction, the RTI, and the heat flux driven buoyancy instability (HBI). In simulations 
with RTI-stable contact discontinuities, the temperature discontinuity spreads due to 
vertical heat conduction. This occurs even for initially horizontal magnetic fields due 
to the initial vertical velocity perturbation and numerical mixing across the interface. 
The HBI slows this temperature diffusion by reorienting initially vertical magnetic 
field lines to a more horizontal geometry. In simulations with RTI-unstable contact 
discontinuities, the dynamics are initially governed by temperature diffusion, but the 
RTI becomes increasingly important at late times. We discuss the possible application 
of these results to supernova remnants, solar prominences, and cold fronts in galaxy 
clusters. 

Key words: instabilities - conduction - MHD - Sun: CME - ISM: supernova rem- 
nants - galaxies: clusters: intracluster medium 



1 INTRODUCTION 

The interface between fluids of different densities can be 
destabilized by gravity or other sources of acceleration. If 
the fluid is accelerated in a direction opposite of Vp, the in- 
terface is unstable to the Rayleigh- Taylor instability (RTI). 
When a heavy fluid is on top of a lighter fluid, the RTI is 
characterized by bubbles of light fluid rising into the heavy 
fluid, and downward spikes of heavy fluid penetrating into 
the light fluid (e.g.. 



produced by AGN in galaxy clusters (e.g., Robinson et al 



2004) 



In some of these contexts, the plasma is hot and dilute. 
If the electron gyroradius is much smaller than the electron 
mean free path, then the heat conduction by electrons is pri- 
marily along the magnetic field |Braginskii 1965). We will 
refer to this effect as anisotropic conduction. Work by, e.g., 



Balbus (20001, Quataert (20081, and McCourt et al. (20111, 



Chandrasekhar 19611. The RTI is po- 



tentially important in a variety of astrophysical situations, 
e.g., in supernova remnants at the contact discontinuity be- 
tween the shocked supernova gas and the shocked interstellar 
mate rial ( |Gull||1973| |Chevalier et al.|[l992] |Jun fe Norman 
19961; in emerging magnetic flux in the solar corona (Isobe 



et al. 2005 Bcr ger et al.|201l[ ); and in buoyant bubbles of gas 



* E-mail:dlecoanet@berkeley.edu 



has shown that convection is strongly affected by anisotropic 
conduction. The normal Schwarzschild condition for connec- 
tive instability, ds/dz < 0, is modified when the conductivity 
time-scale is much shorter than the buoyancy time-scale. In 
this 'fast conduction' limit, plasmas are unstable to a heat 
flux driven buoyancy instability (HBI) when g • VT < or 
are unstable to the magnetothermal instability (MTI) when 
g ■ VT > 0. Although the HBI acts to reduce the effective 
conductivity in the direction parallel to gravity, the MTI 
produces vigorous convec tion (e.g., |Parrish fc Stone||2007| 
Parrish fc Quataert|2008) |McCourt et al.|20TiJ 
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Our goal in this paper is to investigate the magnetized 
RTI when anisotropic conductivity is important. |Balsara| 
et al. ( 2008a I and Balsara et al. ( 2008b I have simulated 



supernova remnants with anisotropic conduction. However, 
their simulations are of an entire supernova remnant, and 
their analysis did not describe how the contact discontinu- 
ity is affected by the combination of the RTI and anisotropic 



conductivity. In fact, in Tilley & Balsara (20061, there does 



not seem to be any RTI type perturbation of the contact 
discontinuity, unlike in, e.g., Chevalier et al. (19921. 



Our analysis will make several simplifying assumptions. 
We will focus our attention only on weak magnetic fields that 



are not dynamically important. Stone & Gardiner (20071 



have shown that stronger magnetic fields prevent fluid mo- 
tions which bend field lines, and that in the strong magnetic 
field limit, the RTI dynamics are dominated by interchange 
motions. Our assumption will be that the fields are strong 
enough to keep the conduction anisotropic, but not so large 
that they are dynamically important. We will also largely 
(though not entirely) neglect the effects of anisotropic vis- 
cosity, which acts to preferentially damp motion along mag- 
netic fields. Previous work has shown that this also leads to 
the RTI being dominated by interchange motions ( |Dong fc| 
Stone|2009| ). 



In addition to studying the RTI, we will also investigate 
the effects of anisotropic conduction on contact discontinu- 



ities that are RTI stable. Birnboim et al. ( 2010 1 propose that 



cold fronts in galaxy clusters can form through the merger of 
shocks. These observed cold fronts are sharp contact discon- 
tinuities, and the intracluster medium (ICM) is sufficiently 
hot and dilute for anisotropic conduction to be important. 
These contact discontinuities are likely to be RTI stable, 



given their small widths. Parrish & Quataert (20081 and 
|Birnboim et al.| ( |2010[ ) suggest that anisotropic conduction 
could help to maintain the small widths of cold fronts in 
clusters, because the HBI acts to suppress heat conduction 
at RTI-stable contact discontinuities. 



The remainder of the paper is organized as follows. 
First, we perform a linear stability analysis including the 
effects of anisotropic conductivity on an RTI-unstable con- 
tact discontinuity with a weak horizontal magnetic field (?|2|. 
Next, we discuss our numerical methods for simulating this 
setup using fully non-linear compressible MHD simulations 
in S|3] In ^ we consider the hydrodynamic problem, i.e., 
without magnetic fields or anisotropic conduction. Section[5] 
studies RTI-stable and RTI-unstable contact discontinuities, 
for a range of initial magnetic field orientations. Finally, in 
S|6]we conclude and discuss the astrophysical implications of 
this work. 



2 LINEAR THEORY 
2.1 Basic equations 

Throughout this paper, we employ the same set of basic 
equations. They are as follows: 



dp 

dt 

dv 

p Tt 

dB 

~dt 



-pV • v, 

-Vp+ —V x (V x 

4-7T 

V x (v x B) , 



V-Q, 



gpz, 



(i) 

(2) 
(3) 
(4) 



where p is the density, v is the velocity, p is the pressure, B 
is the magnetic field, g is the strength of gravity, and z is 
the unit vector in the z direction. T is the temperature, s is 
the entropy, and Q is the heat flux given by 



Q = -bftb • VT, 



(5) 



where b is the unit vector in the B direction, and k is the 
thermal conductivity. We also use d/dt to denote the full 
time derivative, i.e., d/dt = d/dt + v • V. We will employ 
the thermal diffusion coefficient \ — Tk/P in lieu of n. For 
simplicity, we will approximate \ to be constant. Through- 
out this paper, we will take the medium to be an ideal gas 
with ratio of specific heats 7 = 5/3, and use units where 
kbprrip = 1, where kf, is the Boltzmann constant, p is the 
mean molecular weight, and m p is the proton mass. In these 
units the equation of state is P = pT, and T has units of ve- 
locity squared. We will use S to denote perturbed quantities 
below, i.e., Sp is the Eulerian perturbation of the density. 
In the linear problem discussed here, we will assume that 
the magnetic field is sufficiently weak that it is dynamically 
unimportant, and thus the Lorentz force can be dropped 
from the momentum equation. However, the Lorentz force 
term is kept in our non-linear simulations. In several simu- 
lations we also include effects due to anisotropic viscosity in 
the momentum and energy equations (see j |5.4| , but we do 
not include anisotropic viscosity in our linear analysis. 

2.2 Background 

We will now describe the background quantities in our linear 
stability problem and the properties of the perturbations. In 
the RTI problem, two layers of different densities and tem- 
peratures are separated by an interface at z = 0. The pres- 
sure is taken to be continuous at the interface. We will per- 
turb this equilibrium and investigate the linear stability of 
these perturbations. We restrict our attention to magnetic 
fields which are initially horizontal. This allows us to use 
a background state with a density discontinuity (and thus 
a temperature discontinuity) which still satisfies the energy 
equation (eqn.|4|. To avoid any conduction within the upper 
and lower layers, we will assume that the background state 
in each layer is isothermal. This means that the background 
pressure and density vary as exp(—z/H) where H is the 
density scale height. H is related to c s , the adiabatic sound 
speed, by H = ^/(■yg). By using isothermal domains above 
and below the temperature discontinuity, we can restrict our 
attention to instabilities caused by the discontinuity, as op- 
posed to conduction-mediated instabilities within the upper 
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and lower domains (e.g., the MTI or HBI). Although p is no 
longer constant in the upper and lower isothermal layers, we 
will still refer to p + = p\ z =+e and p- = p\ z= - e . 

The pressure is continuous at z = so the pressure 
jump, A [P], is zero, where we define 



a [/] = /!,=+« -/|, = . 



(6) 



for any function /. This implies A [p] = -A [T]. The RTI 
corresponds to A [p] > 0, which yields A [T] < 0, which we 
would expect to also be unstable to the MTI ( |Balbus|2 000). 
However, states with A [p] < 0, and thus A [T] > 0, are 
RTI stable and we would expect to be MTI-stable as well 
(but for non-horizontal fields, potentially HBI-unstable) . We 
summarize the instabilities associated with each tempera- 
ture and density jump, for different initial magnetic field 
geometries, in Fig. [I] 

We assume that all perturbation quantities vary in the 
horizontal and time directions as 



Sp(x, y, z, t) — Sp(z) exp (ik x x + ik y y — iut) . 



(7) 



Although we Fourier decompose the perturbations in the 
horizontal directions, we do not Fourier decompose the 
perturbations in the vertical direction. We will denote by 
k± — ^Jk% + ky the wavenumber perpendicular to gravity 
(not perpendicular to the magnetic field). 



2.3 Boussinesq limit 



Following e.g., |Chand rasckhar (19611 or Shivamoggi (20081, 
we will derive the dispersion relation in two steps. First, we 
use the equations for the perturbations in the upper and 
lower layers to derive the vertical structure of the perturba- 
tions for z^O. Second, we apply the boundary conditions. 
The first three are that the vertical kinetic energy decays 
to zero at z — > ±oo, and that the vertical velocity, Sv z is 
continuous at z = 0. The vertical derivative of the vertical 
velocity, d z 5v z , is in general discontinuous at z = 0. We thus 
take the limit of the equations as z — > to derive a jump 
condition for d z 5v z , which will lead to a dispersion relation. 

First we will simultaneously employ the local and 
Boussinesq limits. The local limit amounts to neglecting any 
derivative of background quantities in the upper and lower 
domains in comparison to the derivative of the perturba- 
tions (Hk± ^> 1). In the Boussinesq limit, 8P/P is small 
in comparison to e.g., 8p/p. In these limits, the momentum 
equation (eqn. [5| gives 



U) V (pSv z ) 



(pSv) = k^giuiSp, 



(8) 



First consider the upper and lower layers, where z / 0. 
One can use either the energy equation or the continuity 
equation to solve for 8p. The energy equation implies 



where 



(icj — cj c ) — + Sv z d z s — 0, 
P 



2 tfl 
5 



(9) 



(10) 



is the conduction time-scale, and we take the magnetic field 
to be in the x direction. The d z s term is proportional to 
1/H, which suggests it might be small in the local limit. It 
can be shown a posteriori that it is consistent to ignore this 



term (using that Hk± 3> 1). When we assume that 8v z d z s 
is small, we get that 8p — in the upper and lower layers. 
Then Equation [8] reduces to 

(-fcl + dl) Sv z = 0. (11) 

This implies 



A+ exp (+k±z) + A- exp (— k±z) 



(12) 



for some amplitudes A+ and A-. This is the same result as 
for the classical RTI analysis (i.e., constant density layers 
with no magnetic fields or conduction), so anisotropic con- 
duction has not affected the vertical structure of the pertur- 
bations. 

The next step is to apply the boundary conditions. The 
boundary conditions at infinity and the continuity of 5v z at 
z = imply 



6v z (z) = 8v z (0) exp (— . 



(13) 



where 8v z (0) is an amplitude. In general we can write 8v z ~ 
exp(k±z), where the plus sign is taken for z > and the 
minus sign is taken for z < 0. The boundary condition at 
infinity requires Re(fc+) < 1/(2H) and Re(fc_) > 1/(2H). 
In this for the classical RTI problem, we have k± = 

Now we will discuss the jump condition on d z 8v z at 
z = 0. Again we must combine the momentum equation, 
Equation|8] with the energy equation. The energy equation, 
in the limit that z — > 0, can be written as 

(iu — u c ) (— icuSp + 8v z d z p) — 0. (14) 

Thus, we have that 

— iuiSp + 5v z d z p — (15) 

regardless of the value of w c - This can be combined with 
Equations [13] and [2] to give the same dispersion relation as 
for the classical RTI problem, 

^ = gk ± (p--p + ) ^ 
P+ + P- 

where A = (p+ — p-)/(p+ + P-) is the Atwood number. 

The key implication of Equation [16] is that in the lo- 
cal, Boussinesq limit the dispersion relation for the RTI is 
independent of the magnitude of the thermal conduction, 
provided that the magnetic field is aligned with the con- 
tact discontinuity. We will now briefly discuss this result. 
In the local limit, d\ogP/dz is much smaller than k±, and 
we always have that A [P] =0. Furthermore, the Boussi- 
nesq approximation implies that SP/P -C Sp/p. These facts 
together mean that the Lagrangian pressure perturbation, 
P' IP is small everywhere, where we use /' to denote the 
Lagrangian perturbation of /. The equation of state implies 
that 



P[ 
P 



! rrit 

p T ' 



(17) 



but since P' / P is small, we have that 

£=-£-=0, (18) 
T p 

where the last equality is due to the continuity equation 
(eqn.[T]). This implies that in the Boussinesq approximation, 
perturbations are isothermal, and thus are unaffected by 
conduction. 
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Figure 1. A schematic of the various instabilities of a contact discontinuity with anisotropic thermal conduction. On the left we show 
the contact discontinuity with a temperature and density of T+ and p_|_ immediately above the interface, and a temperature and density 
of T_ and p_ immediately below the interface. On the right we summarize the instabilities associated with each temperature and density 
jump, for either horizontal or vertical initial magnetic fields. 



The assumption that A [P] = at the interface is a key 
part of the RTI problem. However, the assumptions that 
8P/P is small and that dlogP/dz is much smaller than k± 
result from the Boussinesq and local approximations. 



2.4 Fully compressible linear theory 

We will now relax the Boussinesq and local approximations 
and solve the linear fully compressible RTI problem includ- 
ing anisotropic thermal conduction. The solutions of the 
fully compressible equations are not isothermal, and thus 
are affected by anisotropic heat conduction. 

As above, we must first use the evolution equations in 
the upper and lower layers to solve for the vertical struc- 
ture of the perturbations. After some algebra, one finds the 
following equations for k±: 



fc+ = 



2H+ - \/ m\ 



1 



+ k 2 , -R 



(1-R) 



H+uj 2 ' 



(19) 



1 =2ir^liHi +k -- R 9 w=-^- R ^'^ 

where H± — T±/g is the density scale height in the upper 
and lower layers, and 



R 



(21) 



The quantity R represents the ratio of the response time- 
scale for pressure perturbations to the response time-scale 
for density perturbations. If lo 2> then R — 3/5 and 
perturbations are adiabatic. If ui c 2> w, then R = 1 and 
perturbations are isothermal. 

Note that there are two roots for each of k+ and fc_ in 
Equations 1 19| fe [20| The physical solution has perturbations 
with finite kinetic energy at infinity, i.e., Re(fc+) < 1/(2H+) 
and Re(fc_) > 1/(2H~). Since the two roots of k+ sum to 
1/(2H+) and the two roots of fc_ sum to 1/(2H-), at least 
one root for each of k± will be physical ( Cunningham et al. 

[2oTTb. 



By integrating the momentum equation from z — — e 
to z — +e and dropping terms of order e, we derive the 



following dispersion relation: 



= gk 2 ± A[p] -w A 



R7 



(22) 



where k z \± e = k±. Along with Equations |19| fc |20| this can 
be used to solve for oj as a function of the various parameters 
of the problem - g,k±,oj c , p± and T±. 

Note that for the fully compressible problem, the val- 
ues of k± depend explicitly on to, so we cannot solve for the 
vertical structure of the perturbations before solving for the 
growth rate - instead we must solve for oj and k± simultane- 
ously. Also, because the quantities inside the square roots in 
Equations [19] & [20] are generally complex, we do not know a 
priori if the square root term will have positive or negative 
real part. Our procedure for finding physical modes is as 
follows. We first pick a root for each of k+ and and then 
search for an uj which satisfies the dispersion relation (eqn. 
22j. Then we check that for this co, Re(fc+) < 1/(2H+) and 
Re(fc_) > 1/(2H-). If these conditions are met, then the 
mode is physical. 

To reduce the dimensionality of our parameter space, 
we nondimensionalize lengths and times by setting k± — 2n 
(i.e., we set the perpendicular length scale to one) and g = 1. 
We also take the smaller of p± to be 1. The condition that 
T± — gH± means that setting T± determines the density 
scale height, as well as determining the sound speed through 
c^ ± = ■yT±. Continuity of pressure at z — requires that 
T+p+ = T-p-. With all this in mind, the problem has the 
following degrees of freedom. First, we must pick if the top or 
bottom of the density discontinuity is at the lower density 
of p = 1, as well as the density on the other side of the 
discontinuity (these together are equivalent to picking an 
Atwood number). Next, we pick the scale height on one of 
the sides, which uniquely sets the scale height on the other 
side. Finally, we must pick uj c . 

For the unstable case, we take p_ = 1. For simplic- 
ity, we will take p+ = 3 (A — 1/2), but other values of 
p+ give qualitatively similar results. We now have only two 
remaining parameters: H+ and uj c . In Fig. [2] we plot the 
growth rate, normalized to cjrti = V 'Agk± as a function of 
H+k± for several values of w c /wrti- Without compressible 
effects, the growth rate would be one when normalized to 
wrti- Note that compressibility effects become important as 
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Figure 2. Growth rates, oj, normalized to ujr.ti = VAgkj_ as a 
function of ii+fcx for various values of anisotropic thermal dif- 
fusion, parameterized by oJc/wrti (the ratio of the conduction 
frequency to the classical RTI growth rate). For (j c /o;riTl > 10, 
we are in the isothermal, or fast conduction, limit and the re- 
sults are nearly equivalent to those shown for oJc/urti = 10. As 
-ff+fcx decreases, compressibility becomes more important and 
w/ojrti decreases. At a fixed H+k±, the growth rate increases as 
uj c increases from zero - an asymptotic value is reached in the 
isothermal, or fast conduction, limit when i^c/ajptxi ^ 1. 



i.e., when H+k±_ ~ 1. 
Compressibility decreases the growth rate of the instability, 
as energy is used to compress the fluid. At large cj c , the per- 
turbations are isothermal (i.e., we are in the fast conduction 
limit), whereas at small cj c , the perturbations are adiabatic. 

We will now discuss the differences between the large 
conductivity (isothermal, fast conduction) and small con- 
ductivity (adiabatic) limits. To make the analysis easier we 
will take the limit in which H+k± is small. This corresponds 
to the limit in which the wavelength of the mode is much 
larger than the scale height of the background. It will turn 
out that w 2 / w rti ~ H+k± as H+kx — > 0, which we will 
check a posteriori. First consider the isothermal limit. When 
R = 1, the equations for the vertical structure (eqn. 19 & 20 1 
become much simpler as the 1 



- R terms are zero. Given that 
uj 1 /uJjiTi ~ H+k±, the most important term in the square 
root when H + k± is small is the 1/4H 2 - term. Thus, in this 
limit, k+ « and « 1/H-. 

The dispersion relation then becomes 



— uj 



p-) 



(23) 



This can be easily solved to give 



,2 tt P+ — P- 

-gk±H+ 



lu c /oj > 1, H+k ± < 1. 



( 24 ) 

We see that cj 2 differs from u\ TI by a factor of H + k±{p+ + 
P-)/p- 

Now consider the small conductivity (adiabatic) limit. 
We will still take the limit in which H+k± is small (i.e., the 
long wavelength limit). This problem is more difficult as the 
1 — R terms in Equation [l9]&[20] are now important and u> 



appears explicitly in the formula for k±. We find that 

(26) 



H-k- 



1 l_2^p+ 

2 V 4 5 w 2p_' 



where cj 2 = <jj 2 /(H+k'±) is a normalized growth rate. The 
dispersion relation can then be written as 



-UJ 2 

5 p+ 



H+k. 



= 1 



P+ 



3—p- 

-U) z 



H-k. 



P+ 



(27) 



This is an equation for u 2 with only one free parameter, 
p-/p+. We were unable to find a general analytic solution for 
arbitrary values of p_/p+. However, the equation can easily 
be solved numerically. For P-/p+ = 1/3 or A = 1/2, we have 
that uj 2 w —0.992. Note that for the isothermal calculation, 
when p_/p+ = 1/3 or A — 1/2 we have uj 2 = —2. 

We generally find that the growth rate for isothermal 
perturbations is faster than for adiabatic perturbations. This 
can be explained by describing how isothermal and adia- 
batic perturbations move in an isothermal atmosphere. The 
ratio of densities between a perturbed fluid element and its 
surroundings remains constant if the perturbation and back- 
ground are isothermal. However, for adiabatic perturbations, 
the density of an upward propagating (low density) pertur- 
bation drops less quickly than the density of the background, 
so the perturbed element is less buoyant as it rises. Similarly, 
an adiabatically falling (high density) perturbation increases 
in density less quickly than the background, so the perturbed 
element feels a smaller downwards gravitational force as it 
falls. 

Increasing the anisotropic conductivity from zero to in- 
finity, the value of the R parameter moves from one to 3/5 
through the complex plane. We have not found any oversta- 
bilities in this problem. The behavior at low H+k± is given 
by Equations 25)27 replacing 2/5 by 1 — R and 3/5 by R. 
The growth rate increases smoothly from the adiabatic value 
to the isothermal value as ui c increases. Most of the change 
in uj as uj c varies occurs when uj c is near the adiabatic and 
isothermal growth rates. 

The most prominent feature of the growth rates as 
H+k± becomes small (i.e., for long wavelength modes) is 
that the growth rates also become very small. One possi- 
ble interpretation is that these modes are very compressible 
and most of the potential energy gained by changing vertical 
position is absorbed by compression. A different interpreta- 
tion involves the background density profile. Assume that 
a mode of wavelength k± is 'sensitive' to a distance about 
1/fcx above and below the interface. If H+k± 3> 1, then 
the density is almost constant above and below the inter- 
face. However, if H+k± <g 1, there are many density scale 
heights that can fit within a length ~ l/k±. In this limit, 
the density jump at z — looks like a small perturbation to 
an almost isothermal atmosphere. In this situation the RTI 
is relatively insignificant because the density jump is small 
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compared to the isothermal stratification, so the growth rate 
should also be small. 



3 NUMERICAL METHODS 

We will now discuss the evolution of the RTI in the non- 
linear regime. We will first discuss our numerical methods, 
and then review and clarify previous results on the hydro- 
dynamic RTI. We then run simulations of the RTI problem 
with conduction, with various magnetic field geometries. We 
will primarily be tracking the non-linear evolution by mea- 
suring the height of the highest buoyant bubble as a function 
of time. 

We perform fully non-linear simulations using the 
Athena code ( |Gardiner fe Stone||2"008| |Stone et aLl|2008[ ), 
which solves the compressible evolution equations in con- 
servative form using a Godunov method. We implement 
anisotropic conductivity as described in |Parrish fc Stone] 
(20051) andlSharma & Hammett] p007l. 



Table 1. Parameters for simulations without conductivity. 



The simulation domain is Cartesian, with the computa- 
tional domain x,y £ [— L/2, +L/2], z £ [— L, +L] for some 
box length L. We normalize our lengths to L (i.e., set L = 1) 
in all our simulations. Simulations are run at three res- 
olutions: 64x64x128, 128x128x256, and 256x256x512. The 
boundary conditions are periodic in the horizontal direc- 
tions. At the top and bottom boundaries, the pressure is 
extrapolated to maintain hydrostatic equilibrium, and re- 
flecting boundary conditions are applied to the remaining 
variables. 

Our initial conditions for temperature are T = T+ for 
z > and T — T- for z < 0. The pressure and density, 
within the isothermal domains, are given by 



P ± (z),p±(s)~exp(-gz/T ± ), 



(28) 



so the scale height within the domains is H± — T±/g. We 
will refer to p± = p\± e . In all of our simulations we take the 
smaller of p± to be one and the larger to be three. The ratio 
H±/L relates the sound crossing time of the box to the RTI 
growth time for perturbations with wavelengths comparable 
to the size of the box. In all of our simulations, we pick 
H/L = 100 and g = 1, where H is the scale height on the 
p = 1 side of the interface. We take 7 = 5/3, so c 2 — 7T = 
500/3 on the side of the interface with p = 1. The sound 
crossing time is w 13, in comparison to the RTI growth time 
for modes the size of the box, w 1.7. Thus, our simulations 
are fairly incompressible. 

The uniform magnetic field is taken to be horizontal, 
vertical, or forming a 45 degree angle with the horizontal. 
The magnetic field strength is set to B/\ // 4tt — 0.0001. For 
the linear RTI, such a field would produce noticeable dy- 
namical effects only on length scales comparable to 5 x 10~ 7 



(Stone & Gardiner 2007 hereafter SG07), which cannot be 
resolved by our simulations. Furthermore, the Alfven ve- 
locity is about 10 _5 c s . Thus, the magnetic field plays no 
dynamical role. 

The conductivity, x, is normalized to yf gL 3 . That is, if 
we take ujrti = y/2TvAg/L and lj c — §x4tt 2 / L 2 , where L is 
the length of the box, then uj c /ujrti = 2(27r) 1,5 /(5V^4)x « 
9x for our problem. We will primarily pick values of x 
smaller than one to understand the interaction of conduc- 
tivity with the RTI. 



Name 


Resolution 


Constant p or T 


SG 


256x256x512 


P 


RTLR 


64x64x128 


T 


RT 


128x128x256 


r 


RTHR 


256x256x512 


T 



The simulations have either constant p layers, as in SG07, or con- 
stant T layers. The simulation with constant p layers is run using 
the same parameters as the hydrodynamic run in SG07. These 
parameters are described at the beginning of j|4 The simulations 
with constant T layers are run as described in f 3] 



We start the simulations with a small perturbation 
to the vertical velocity. We set the velocity to AqR(1 + 
cos(nz/L)), where Ao the amplitude and R is a random 
number between —1 and +1, as in SG07. The z dependence 
is chosen so there is no vertical perturbation at the top and 
bottom boundaries. We take Ao = 0.05, which implies that 
the maximal velocity perturbations (of size 2Ao = 0.1) are 
~ 0.01c s . The initial perturbations are large enough to start 
the simulations close to the non-linear regime. In Appendix 
A we discuss simulations with smaller initial perturbations 
which probe the linear regime discussed in S|2] 



4 RTI WITHOUT CONDUCTIVITY 

Before presenting the results of our simulations, we briefly 
discuss previous results on the non-linear evolution of the 
RTI as presented in Dimonte et al. (2004 1. Non-linear effects 



become important for perturbations of wavelength A = 2-7r/fc 
when the perturbation has moved the interface a height ~ A 
up or down. The ascending light fluid is described as a bub- 
ble, whereas the descending heavy fluid is described as a 
spike. The main diagnostic we will use is the height of the 
highest bubble as a function of time. The results for the 
depth of the lowest spike as a function of time are qual- 
itatively similar. In laboratory experiments and numerical 
simulations, the height h is reported to grow as 



h = othAgt , 



(29) 



where A is the Atwood number and aj, «0. 04-0. 08 is a di- 
mensionless number arrived at experimentally. Numerical 
simulations produce values of oj, around a factor of two 
lower than physical experiments - this is thought to be due 
to increased diffusion at the interface in numerical simula- 
tions ( |Dimonte et al.|2004 i. 

To provide some context for the calculations with 
anisotropic conduction to follow, we now make a detailed 
comparison with the hydrodynamic RTI simulation of SG07. 
We have run a numerical simulation using the exact same 
parameters as reported in SG07^Ve compare this to a sim- 
ulation using the initial conditions described in ^3] that we 
will use for our simulations with anisotropic thermal con- 
duction. 



1 Specifically, we ran the problem with constant density layers 
with P = 3/5 at z = at a resolution of 256x256x512, without 
magnetic fields or conduction, and took L = g = 0.1. There is 
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In previous work on the RTI with constant density lay- 
ers, the height of the highest bubble is typically defined us- 
ing the fraction of heavy fluid (fh) and the fraction of light 
fluid (fi = 1 — fh) within a volume ( Dimonte et al.|[2004 |. 
If we use p+/p- = 3, then p = 2fh + 1. At t = 0, fh = 1 
above the temperature discontinuity. Thus, many define the 
height of the highest bubble to be the highest point at which 
the horizontally averaged fh deviates from one by a small 
amount. We will denote horizontal averages of quantities us- 
ing (■), e.g., the horizontal average of fh is (fh)- In SG07, 
the height is defined to be the highest point at which (fh) 
is less than 0.985. This corresponds to the highest point at 
which the density is less than p — 2.97. Using a lower cut- 
off density (or equivalently, a smaller value of fh) does not 
qualitatively change our results. 

In Table[T]we list our simulations which have no conduc- 
tivity. In run SG we use the parameters reported in SG07. 
In the simulations beginning with RT, we use the initial con- 
ditions described in S|3] We use the suffix LR to denote low 
resolution runs with resolutions of 64x64x128 and HR to de- 
note high resolution runs with resolutions of 256x256x512. 
RT is our fiducial model which we will compare to our sim- 
ulations with anisotropic conduction. 

In the runs labelled RT and the simulations with 
anisotropic conduction presented in S|5] the density is a func- 
tion of height, but the temperature is initially constant in 
the upper and lower layers. Thus, we will define the height 
to be the highest point at which (T) differs from T + by more 
than 1 per cent. This is equivalent to the definition of height 
used by SG07 provided that the pressure does not change 
significantly through the simulation. In run SG, we use the 
same definition of height as SG07. 

We plot the height for runs SG (using the same pa- 
rameters as SG07) and RTHR (our high resolution simu- 
lation without conduction) in the top panel of Fig. [3] The 
two simulations give similar results. Interestingly, they both 
differ somewhat from the data in SG07 (their fig. 8). It is 
unexpected that run SG has a different height than the cor- 
responding simulation in SG07, as we have used the exact 
same parameters and a newer version of the same code as 
reported in SG07. Nevertheless, our results are qualitatively 
similar. The differences between our results and those of 
SG07 are presumably due to a small difference in initial 
conditions or parameters of the run that we have not been 
able to identify. 

In the upper panel of Fig. [3] we plot h/L as a function 



of Agt 2 / L, as is often done in the literature (e.g., |Dimonte 



et al. 12004 1. We expect to have h - 
straight line in this figure. Near t 



a.hAgt , so we expect a 



0, the slope is changing 
rapidly, but the lines are fairly linear at later times. However, 
it is very difficult to tell from a plot of h vs. t 2 whether or 
not the height is truly increasing as t 2 . 



some ambiguity regarding the initial perturbation in SG07. They 
state that their initial perturbation is smoothed toward the verti- 
cal boundaries, but then write that v z (z) = Aoii(l-|-cos(27rz/L)), 
where Aq is an amplitude and R is a random number between — 1 
and +1. Note that v z is maximal at the boundaries if the vertical 
boundaries are at z = ±L. We believe this is a typo and they 
actually set v z (z) = AqR(1 + cos(-7r,z/L)), so that v z (±L) = 
— we use this functional form for the perturbation. We also take 
A = 0.005. 




CD 
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t/(L/Ag) 1/2 
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Figure 3. The height of the highest bubble (defined in jj4]l as a 
function of time in RTI simulations without conduction. In the 
top panel the height is plotted against t 2 . In the bottom panel the 
height is plotted on a log-log plot against t. We show the results 
from our high resolution simulation without conduction (solid 
line; RTHR in Table [TJ and our simulation using the parameters 
of SG07 (dashed line; SG). The two dotted lines show power law 
relations with exponents of one and two. Although the height vs. 
t 2 line in the upper panel might look fairly linear, the bottom 
panel shows that the height grows linearly with t at later times. 



To test whether h is in fact increasing as t 2 , we plot 
vs. i on a log-log plot in the bottom panel of Fig. [3] The 
dotted lines show power laws with exponents of one and two. 
Although the two simulations are initially roughly consistent 
with h ~ t 2 , for times past t/ \/L/Ag « 1.3, the results 
are more consistent with the height increasing linearly with 
time. The numerical results from SG07 are also consistent 
with the height increasing linearly with time. 

In Fig. [4] we plot, for three different resolutions, the 
height normalized to L as a function of Agt 2 normalized to 
L. There is no significant agreement between the simulations 
at different resolutions. However, in the bottom panel of Fig. 
[4] we show the height normalized to Ddom as a function of 
Agt 2 normalized to Ddom, where Ddom is the size of the first 
dominant bubbles and spikes in each simulation. There is 
much closer agreement between the runs at different reso- 
lutions using this normalization. This effect was previously 



noted by, e.g., Dimonte et al. (20041 



The statement that h should increase as t 2 is often 
motivated by dimensional analysis: if the only dimensional 
parameter in the problem is g, then we must have that 
h = Cgt 2 for some dimensionless coefficient C (which is 
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t 2 /(L/Ag) 

400 E 
E 300 - 




200 400 600 800 1000 1200 
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Figure 4. The height of the highest bubble (defined in ij4j| as 
a function of time for our RTI simulations without conduction 
at different resolutions. In the upper panel the height and Agt 2 
are normalized by the length of the box L. In the lower panel, 
the height and Agt 2 are normalized by Ddom, the size of the 
first dominant bubbles. We plot our results for our low resolution 
run (solid line; RTLR in Table [TJ, our fudicual run (dashed line; 
RT), and our high resolution run (dot-dashed line; RTHR). We 
estimate D^om to be L/32, or about 8 cell widths for our high 
resolution simulation. We then take -Ddom = L/16 and L/8 for 
our fiducial simulation and low resolution simulation respectively. 
The results line up well when normalizing lengths to Dd om , but 
there is little agreement between the different resolutions when 
normalizing lengths to L. This indicates that the simulations do 
not notice the size of the box, but instead remember the size of the 
first dominant modes. The flat portion of the curves correspond 
to when the perturbations have hit the top of the box in each of 
the simulations. 



related to the Atwood number). However, our simulations 
show there is another dimensional number in the problem: 
Ddom- When we increase the resolution of our simulations by 
a factor of two in each direction, Ddom decreases by a factor 
of about two. The bottom panel of Fig. [4] shows that the 
relevant length scale for the non-linear evolution is Ddom, 
and indicates that the simulations remember the size of the 
initial dominant modes very late into the non-linear regime. 

With another dimensional number in the problem, 
Ddom, we can make a new dimensionless parameter, 
Ddom/(gt 2 ). With this additional dimensionless parameter, 
the height is no longer required to increase as gt 2 . Instead, 
for instance, the height could increase as t^/D^^g. Physi- 
cally, this would correspond to a bubble of fixed size rising at 




1.0 



1.0 



Figure 5. Isosurfaces of temperature at T = 95 and T = 38 for 
the high resolution simulation without any explicit conductivity 
(RTHR in Table [TJ. The upper (lower) layer is initialized with a 
temperature of 33 (100). The upper panel is at time t/ ^/L/Ag Ri 
1.8, and the lower panel is at time t/^L/Ag £3 4.2. Temperature 
contours are also shown on the faces of the simulation domain. 
The height of the highest bubble is increasing linearly at both 
times shown (see Fig. [3j, consistent with a bubble of fixed size 
rising at its terminal velocity. However, the size of the bubbles 
is increasing, which would yield larger terminal velocities at later 
times. 

its terminal velocity. This functional form fits our numerical 
results, although the size of the bubbles in our simulation 
increase as a function of time (see Fig. [5| . 

Although we do not have a satisfactory explanation for 
why the height increases linearly with time, it is clear that 
our results are not consistent with the commonly repeated 
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result that the height increases quadratically with time. We 
have also shown that the results of SG07 are also not con- 
sistent with the height increasing quadratically with time 
(although SG07 did not make this point explicitly). This 
shows that care must be taken when determining the func- 
tional form of h(t). 

With this in mind, we will now include anisotropic con- 
ductivity in our calculations. We will be making frequent 
comparison with our fiducial RTI simulation without ex- 
plicit conduction (RT in Table [T]), as well as with the results 
shown in Figures [5Jf5] 



5 RESULTS WITH ANISOTROPIC 
CONDUCTION 

Anisotropic conduction enables diffusive spreading of the ini- 
tial temperature discontinuity. Even horizontal fields induce 
diffusive spreading - the small v z from the initial pertur- 
bation produces a small vertical field component. This pro- 
duces a small amount of perpendicular temperature diffusion 
when the interface between the two constant temperature 
layers is within a cell. To distinguish mixing of temperature 
from mixing of material on either side of the interface, we 
will introduce two new definitions of height. We will call the 
definition used in fj4] the height of the temperature mixing 
layer, or hr- As before, it is defined to be the highest point 
at which (T) differs from T+ by more than 1 per cent. We 
will also sometimes calculate the depth of the temperature 
mixing layer, or dx- Similarly, cLt is the lowest point at which 
(T) differs from T_ by more than 1 per cent. 

In addition, we initialize the simulations with two pas- 
sive scalar fields which we will call Cl and Cu- The sim- 
ulations start with Cl = 1, Cu = for z < and with 
Cl = 0, Cu = 1 for z > 0. We define the two heights 
of the composition mixing layer, hc L and hc u , to be the 
highest points for which (Cl) or (Cu) differ from their ini- 
tial value by more than 0.015. That is, hc L is the highest 
point at which (Cl) > 0.015 and hc v is the highest point 
at which (Cu) < 0.985. When there is no temperature diffu- 
sion, h,T = hc L = hcu provided that p+ = 3 and p„ = 1, so 
the height of the temperature mixing layer, the two heights 
of the composition mixing layer, and the height of the high- 
est RTI bubble all coincide. 

In Tabled we list the simulations with conduction that 
we will discuss in this paper. We studied interfaces with 
RTI-stable or RTI-unstable jumps in temperature, as well 
as magnetic fields which are initially horizontal, vertical, or 
at a 45 degree angle to the horizontal (see Fig.[T|. The names 
of simulations with initially horizontal magnetic fields start 
with 'H' in Table |2j whereas the names of simulations with 
initially vertical magnetic fields start with 'V.' The two sim- 
ulations with an initial magnetic field at a 45 degree angle 
to the horizontal have names starting with 'A.' Simulations 
which are RTI stable have an 'S' as the second character in 
their name. We vary the thermal diffusivity, which is given 
by X = 10~ n , where n is the number in the name of the 
simulation (in simulations with initially skewed magnetic 
fields, x = 2 x 10 _n ). A diffusivity of 10~™ corresponds 
to Wc/wrti = 9 x 10~", where uj c is conduction frequency 
over the simulation domain, and cjrti is the RTI growth 
rate for modes the scale of the simulation domain. Some 



Table 2. Parameters for simulations with conductivity. 



Name Resolution lu c /uj^i B dir. RTI stability 



HO 


128x128x256 


9 


hor. 


unstable 


HI 


128x128x256 


0.9 


hor. 


unstable 


H2 


128x128x256 


0.09 


hor. 


unstable 


H2HR 


256x256x512 


0.09 


hor. 


unstable 


HS0 


128x128x256 


9 


hor. 


stable 


HS1 


128x128x256 


0.9 


hor. 


stable 


HS2 


128x128x256 


0.09 


hor. 


stable 


HS2HR 


256x256x512 


0.09 


hor. 


stable 


VI 


128x128x256 


0.9 


ver. 


unstable 


V2 


128x128x256 


0.09 


ver. 


unstable 


V3 


128x128x256 


0.009 


ver. 


unstable 


V3HR 


256x256x512 


0.009 


ver. 


unstable 


VS2 


64x64x128 


0.09 


ver. 


stable 


VS3 


64x64x128 


0.009 


ver. 


stable 


VS3HR 


128x128x256 


0.009 


ver. 


stable 


VS4 


64x64x128 


0.0009 


ver. 


stable 


VS1I 


128x128x256 


0.9 (i) 


ver. 


stable 


VS2I 


128x128x256 


0.09 (i) 


ver. 


stable 


VS3I 


64x64x128 


0.009 (i) 


ver. 


stable 


VS3IHR 


128x128x256 


0.009 (i) 


ver. 


stable 


VS4I 


128x128x256 


0.0009 (i) 


ver. 


stable 


A3 


128x128x256 


0.018 


45° 


unstable 


AS 2 


128x128x256 


0.18 


45° 


stable 



The simulations are run as described in |3] The conductivity is 
anisotropic except for the simulations labelled (i). uj c is the con- 
duction time-scale across the simulation domain and ojrti is the 
RTI growth rate for modes with wavelengths comparable to the 
size of the simulation domain. If the simulation is listed as RTI- 
unstable, then p+ = 3 and p_ = 1. If the simulation is listed as 
RTI-stable then p+ = 1 and p— = 3. 

of our simulations employ isotropic instead of anisotropic 
conduction - these have the letter T in their names. Most 
of our simulations are at a resolution of 128x128x256 (the 
simulations with RTI-stable temperature jumps and vertical 
magnetic field are run for longer, so they use a resolution of 
64x64x128). However, to test resolution effects, we ran some 
simulations at a resolution of 256x256x512 (we ran the sim- 
ulations with RTI-stable temperature jumps with vertical 
magnetic fields at a resolution of 128x128x256) - these high 
resolution runs are denoted by 'HR.' 

5.1 Horizontal magnetic field 

Systems with horizontal magnetic fields and anisotropic con- 
duction are sometimes susceptible to the MTI. The MTI and 
RTI both occur when p+ > p~. We will start by discussing 
the case where p+ < p_. There is neither MTI nor RTI in 
this case, so the dynamics are dominated by diffusion per- 
pendicular to the magnetic field lines due to a combination 
of small perturbations to the field orientation and numerical 
mixing across the interface. If a simulation is run with no 
initial perturbation, we find that the code holds the equilib- 
rium exactly. 

5.1.1 RTI-stable and MTI-stable (p+ < p~) 

Because there is neither MTI nor RTI the distinction be- 
tween the height of the temperature mixing layer (/it) and 
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the composition mixing layer (hc L and hc v ) is very impor- 
tant. In Fig . |6| we plot (T), (Cl), and (Cu) as a function 
of z at t/WL/Ag ~ 3.5 for our high resolution run with 
Wc/wrti = 0.09 (HS2HR in Table [5}. Several effects can be 
seen in this figure. First, the vertical position of the inter- 
face (about where (Cl) and (Cu) intersect) is above z = 0. 
Because the interface is moved upwards, the upper layer 
is compressed, causing Cu to build up above the interface. 
Note that there is very little movement of material across 
the interface - the gradient of (Cu) remains large through- 
out the simulation. Similarly, (Cl) does not spread upwards 
across the interface. However, because the interface is mov- 
ing upwards, Cl is spreading out below the interface. In 
fact, the depth associated with (Cl) is equal to the depth 
associated with (T). Also notice that our metrics for the 
height of the composition mixing layer are sensitive to both 
the width of the composition mixing layer and the vertical 
position of the centre of the mixing layer. Thus, the size of 
the mixing layer is actually being overestimated by hc L and 
hc v , despite remaining small throughout the simulations. 
The results depicted in Fig. [6] are similar for other times - 
at earlier times, the interface is closer to z — 0, whereas at 
later times, the interface is at a higher vertical position. 

The origin of the upwards motion of the interface can 
be understood as follows. Because of the small initial per- 
turbation to v z and numerical mixing across the interface, 
there is some vertical conduction of heat. The perpendicu- 
lar conduction is proportional to the amplitude of the initial 
perturbation which adds a small vertical component to the 
magnetic field. This vertical heat conduction will produce a 
vertical pressure force in the following manner. Consider the 
dynamics at the beginning of the simulation and very close 
to the interface at z = 0. Initially the temperature above 
the interface is higher than the temperature below the inter- 
face, while the pressure is about equal above and below the 
interface. Dynamically nothing is happening because the in- 
terface is stable to dynamical perturbations. However, there 
is some vertical heat conduction, so on a temperature diffu- 
sion time-scale, heat is transferred from above the interface 
to below the interface. Since the temperature above the in- 
terface drops, the pressure above the interface also drops. 
Similarly, the pressure below the interface increases. This 
pressure difference pushes the interface upwards. 

In Fig. [jj we plot the height as a function of time for 
several different values of Wc/wrti- hc L and hc v are only 
plotted for the high resolution simulation (HS2HR in Ta- 
ble [2| , but these heights are similar for the other simula- 
tions. As the conductivity increases, all the height metrics 
increase, although ha L and hc v remain about an order of 
magnitude smaller than hr- Note that hr is slightly larger 
at higher resolutions. This might be because the upwards 
pressure force due to conduction is larger for smaller cell 
widths, as the initial gradient is sharper. We have defined 
hc L and hc v to be the points at which (Cl) and (Cu) dif- 
fer by 1.5%, which corresponds to a 1 per cent temperature 
variation if p+ = 3. However, in these simulations p+ = 1, 
so if T, Cl, and Cu were all mixed by the same amount, we 
would expect hc L = hc n to be larger than hr, instead of 
much smaller than hr as shown in Fig. [jj This highlights 
the fact that the dominant effect in these RTI-stable sim- 
ulations is simple temperature diffusion which increases hr 
much more than hc L or hc v ■ 



A 

h- 

V 



20 
00 
80 
60 
40 





ji 


\ 


i y" 

i / 

1/ 






<T> 








<C L > 




h 


< Cu > 










ji 



1 .5 



o 
V 

a" 

0.5 



-1.0 -0.5 



0.0 



0.5 



0.0 




Figure 6. Profiles of (T), (Cl) and (Cu) for our high resolution 
RTI-st able si mulation with £J c /ojr T i = 0.09 (HS2HR in Tabled 
at t/\/ L/Ag sa 3.5. Cl (Cu) is a passive scalar initialized to 
one (zero) for z < and zero (one) for z > which track the 
mixing of density, rather than the mixing of temperature (which 
is affected by heat diffusion). The vertical dotted line is at z = 0. 
The interface (about where (Cl) = (Cu)) has moved upwards 
from 2 = 0, compressing Cu above the interface and spreading 
out Cl below the interface. This upwards motion is due to a 
pressure force caused by the small vertical heat diffusion in the 
simulation. 
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Figure 7. Heights (defined in jj5jl as a function of time for vari- 
ous simulations with horizontal magnetic fields and an RTI-stable 
temperature jump. The solid and dotted lines are plots of 
the height of the temperature mixing layer. The highest solid line 
correspond to the run with Wc/ojrti = 9 (HS0 in Table [2), the 
middle solid line correspond to the run with with w c /ojrti = 0.9 
(HS1), and the lowest solid line and dotted line correspond to 
the two runs with u) c /u)-^ri = 0.09 (solid line is at the fiducial 
resolution, HS2; dotted line is the high resolution run, HS2HR). 
The dashed and dot-dashed lines show hc L and hc u respectively 
for simulation HS2HR; these are two metrics for the height of the 
composition mixing layer. If there was no perpendicular temper- 
ature diffusion, all heights would stay at zero for the entire simu- 
lation. Although little material diffuses across the interface, hc L 
and hc u increase slowly due to a pressure force which raises the 
interface. 



We also tracked the amplification of magnetic energy 
in the simulations. The results are shown in Fig. [8] In this 
figure, we plot the growth of magnetic energy for all our 
high resolution simulations. For comparison, we include the 
high resolution simulation without any explicit conductivity 
(solid line). In the RTI- unstable simulation without any ex- 
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Figure 8. The magnetic energy, |B| 2 /87r, normalized to the ini- 
tial magnetic energy as a function of time. We show results from 
our high resolution simulation without any explicit conductivity 
(solid line; RTHR in Table [l}; our high resolution simulations 
with initially horizontal magnetic fields with ui c /u>-[iTl = 0.09, 
one which is RTI-stable (dotted line; HS2HR in Table [2), the 
other which is RTI-unstablc (dashed line; H2HR); and our high 
resolution simulations with initially vertical magnetic fields with 
^c/^rti = 0.009, one which is RTI-stable (triple-dot-dashed line; 
VS2HR), and one which is RTI-unstable (dot-dashed line; V2HR). 
The HBI-unstable simulation, VS2HR, is run for much longer 
than the other simulations; for convenience, we rescale the time 
variable using a scaling factor S which is equal to 12 for simulation 
VS2HR and 1 for all other simulations. The RTI-unstable simu- 
lations exhibit exponentially growing magnetic energies, whereas 
the RTI-stable simulations show linearly increasing magnetic en- 
ergies. 

plicit conductivity, the magnetic energy grows exponentially 
for most of the simulation. This is because the magnetic field 
remains dynamically unimportant for the entire simulation. 
In our high resolution RTI-stable simulation with a horizon- 
tal initial magnetic field, there is some initial linear growth 
in the magnetic energy which saturates by tj \J Lj Ag « 2; 
the magnetic energy increases by a factor of four. This is 
likely due to the initial velocity perturbation. Unsurpris- 
ingly, the growth is significantly smaller than in the RTI- 
unstable simulation, although the two cases have the same 
initial magnetic field. 

These results highlight the essential role of perpendic- 
ular thermal conduction in our RTI-stable simulations with 
horizontal magnetic fields. Without the perpendicular con- 
duction, there would be no vertical temperature diffusion, 
so hr would stay close to zero for the entire simulation. The 
interface would also not rise so hc L and hc v would also stay 
close to zero. In any astrophysical context, however, there 
will always be some diffusion across the interface. Temper- 
ature diffusion across the interface could occur due to some 
small perpendicular diffusion, or because the magnetic field 
is not completely parallel to the interface. In these cases 
the system will evolve as shown in Figures [6] & [7] although 
the time scale for diffusion will depend on the source and 
strength of diffusion across the interface. 

5.1.2 RTI-unstable and MTI-unstable (p+ > p-) 

When the density above the interface is larger than the den- 
sity below it, the contact discontinuity is RTI-unstable. We 
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Figure 9. Heights and depths of the temperature mixing layer 
(defined in as a function of time for various i^c/ojrti for sim- 
ulations with horizontal initial magnetic fields. The solid lines 
show hx for our simulations which are RTI-unstable, whereas the 
dotted lines show — dj- for our simulations which are RTI-stable. 
The top pair of blue lines, which are almost overlapping are for 
Wc/ojrti = 9 (HO and HS0 in Table[2]|. Below these, the red pair 
of lines, which are also almost overlapping are for lu c /u)rti = 0.9 
(HI and HS1). The lowest pair of lines, which are in yellow, cor- 
respond to cj c /to>RTl = 0.09 (H2 and HS2). The bottom solid 
line shows hx for our fiducial simulation without explicit conduc- 
tivity (RT in Table [l}. The height for the RTI-unstable runs is 
about equal to the depth for the RTI-stable runs when conduction 
dominates over the RTI, which is for most of the simulation when 
lJc/cjrti = 9 and 0.9. However, for u> c /u)-p,Ti = 0.09, the RTI be- 
comes important around t/^/L/Ag Ri 1-2, and the height in the 
RTI-unstable simulation increases in a similar way to the height in 
the simulation without conductivity. The results are qualitatively 
similar at higher resolution. 

plot hr for the RTI-unstable simulations for several values 
of luc/oj-rti as solid lines in Fig. [9] We also plot the depth of 
the temperature mixing layer, dr, for the RTI-stable simu- 
lations for several values of Wc/wrti as dashed lines. Using 
these two sets of lines, we can determine when the increase 
in height is due to vertical heat conduction and when it is 
due to RTI non-linear mixing. Both the depth of the tem- 
perature mixing layer in the RTI-stable case and the height 
of the temperature mixing layer in the RTI-unstable case 
describe how temperature penetrates into the layer with 
T = 100/3. If vertical heat conduction is more important 
than the RTI, then we would expect hr for the RTI-unstable 
case to equal —dr for the RTI-stable case. Indeed, we see 
that hr is very close to — for the entire simulation when 
oj c /lokti = 9 and 0.9. This indicates that vertical heat con- 
duction dominates the RTI for these diffusivities. However, 
when cj c /ujFrri = 0.09, the vertical heat conduction only 
dominates until tj yj L/Ag w 1 — 2. After this time, the RTI 
becomes more important, and the height increases as in our 
fiducial simulation with no explicit conduction. 

In general, vertical heat conduction is always faster than 
the RTI early in the simulation because the conduction time 
over the length scale associated with the discontinuity is 
nearly zero. Thus, the temperature is essentially unaffected 
by the RTI at early times. However, the RTI does continue 
to grow, even in the diffusion-dominated regime early in the 
simulations. The RTI bubbles grow faster than heat dif- 
fuses, so the RTI always wins at late times provided the 
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Figure 10. Heights (defined in ij5]l as a function of time for our 
high resolution RTI-unstable run with an initially horizontal mag- 
netic field, and with cj c /u HTI = 0.09 (HS2HR in Table[2). We plot 
the height of the temperature mixing layer, hx (solid line), as well 
as the two heights of the composition mixing layer, h.Q L and hc v 
(dashed and dot-dashed lines respectively). When conduction is 
more important than the RTI, h@ L stays small, although hc v 
is about equal to hx- This is because the profiles of (T), (Cl), 
and {Cjj) look similar to those in Fig. [6] except with z — > —z 
and Cl <-> Cjj ■ When the RTI becomes important at around 
t/y/L/Ag Si 1 — 2, everything becomes well mixed, and all three 
height metrics are about equal. 



simulation domain is large enough. The transition to the 
RTI-dominated regime occurs when the height of the RTI 
bubbles reaches the height of the temperature mixing layer 
due to conduction alone. If the RTI was occurring indepen- 
dently of the conduction, this would occur when hr in our 
fiducial simulation equals —dr for the RTI-stable case. This 
is broadly in agreement with Fig. [9] 

However, the height after the transition to the RTI- 
dominated regime is not equal to the height in our fiducial 
simulation which has no explicit conductivity. The simula- 
tions with anisotropic conduction are more dissipative so the 
height increases more slowly than when there is no explicit 
conductivity. This is the principal effect of anisotropic con- 
duction for the horizontal field, RTI-unstable case, in the 
RTI-dominated regime. The increase in dissipation is due 
to magnetic field lines penetrating the bubbles and diffus- 
ing their temperature into the surrounding medium. This 
decreases the buoyancy of the bubbles and slows their rise. 

In the RTI-stable case (p+ < p~), the height of the 
composition mixing layer and the height of the temperature 
mixing layer were fairly decoupled (see Fig. [7}. We plot our 
three height metrics in Fig. |10| fo r our high resolution RTI- 
unstable run (RT2HR in Table |2| . We see that initially hc v 
is about equal to hr, whereas hc L is much smaller. At later 
times, once the RTI becomes important, all three height 
metrics become about equal. 

At early times, when conduction is more important than 
the RTI, the evolution is very similar to the RTI-stable case. 
However, because the high temperature layer is on the bot- 
tom rather on top, the heights in the RTI-unstable case 
match the depths in the RTI-stable case, and Cl and Cu 
are switched. The profiles of (T), (Cl), and (Cu) look very 
similar to those shown for the RTI-stable case in Fig. [6] un- 
der this identification. As can be seen in Fig. [6j the depth 
associated with (Cl) (corresponding to hc u for the RTI- 



unstable simulation) is very close to dr (corresponding to 
fir for the RTI-unstable simulation). However, the depth 
associated with (Cu) (corresponding to hc L for the RTI- 
unstable simulations) is small in comparison to dr- 

In Fig. [TT] we plot isosurfaces of temperature for the 
RTI-unstable simulation with oj c /wrti = 0.09 at two times. 
The upper panel shows the temperature at early times when 
the RTI is just about to start setting the height, whereas the 
lower panel shows the temperature at late times when the 
RTI is more important than conduction. The presence of 
some vertical conduction is clear when comparing the upper 
panel of Fig.[TT]to the upper panel of Fig. [5] which plots iso- 
surfaces of temperature at the same time for a run without 
any explicit conduction. The small scale perturbations on 
the surface of the T = 38 surface are due to a combination 
of uneven vertical conduction - as the vertical conduction 
depends on the size of the initial perturbation - as well as the 
tops of RTI bubbles which are growing despite the diffusion. 
The bubbles eventually grow faster than the temperature 
mixing layer diffuses, and the dynamics become dominated 
by the RTI (as seen in the lower panel). However, conduc- 
tion is still important at late times, as the bubbles do not 
rise as quickly as in simulations with no explicit conduction 
(see Fig. [9}. This can also be seen in the lower panel of Fig. 
[Til - the RTI bubbles are much wider than the RTI bubbles 
in our simulation with no explicit conductivity (Fig. [HJ. 

We find that the presence of the MTI does not signifi- 
cantly alter the RTI dynamics. If we allow our simulations 
to run long enough that the RTI bubbles hit the top and 
bottom walls, we begin to see convection. If the conduction 
time is shorter than the buoyancy time, then this convection 
is effectively driven by the MTI ( |McCourt et al.|2011[ ). How- 
ever, even without conduction, the simulations are convec- 
tively unstable at late times, so the MTI does not introduce 
qualitatively different dynamics into the problem. 

There is significant growth in magnetic energy for the 
RTI-unstable simulations (see Fig. |8|. Just as the magnetic 
energy grows exponentially in the RTI-unstable simulation 
without any explicit conductivity (solid line), the magnetic 
energy also grows exponentially in the RTI-unstable sim- 
ulation with an initially horizontal magnetic field (dashed 
line). The growth rates are similar for both cases, although 
the growth in the simulation with anisotropic conduction is 
somewhat slower (i.e., the exponential growth starts later). 
This is likely because bubbles rise more slowly in the sim- 
ulations with anisotropic conductivity. By the end of both 
simulations, the magnetic energy has amplified by about 10 3 . 

5.2 Vertical magnetic field 

We will now describe simulations initialized with a verti- 
cal magnetic field. Simulations with p+ < p_ are still RTI- 
stable, but are now HBI-unstable, yielding new dynamics. 
The simulations with p+ > p_ are HBI-stable and RTI- 
unstable and are similar to the RTI-unstable simulations 
with horizontal fields. These simulations are not initially in 
equilibrium, as there is a vertical temperature discontinu- 
ity and vertical heat conduction. Physically, the simulations 
should be viewed as probing length scales much longer than 
the initial width of the spread in temperature. On these long 
length scales, the temperature is initially very close to dis- 
continuous. 



© 0000 RAS, MNRAS 000, 000-000 



Ray leigh- Taylor Stable and Unstable Contact Discontinuities with Anisotropic Conduction 13 

1 .0 




1.0 



1.0 



Figure 11. Isosurfaces of temperature at T = 95 and T = 38 
for the high resolution RTI-unstable simulation with an initially 
horizontal magnetic field and with cj c /cjrti = 0.09 (H2HR in 
Table pjl. The upper panel is at time t/y/L/Ag ~ 1.8, and the 
lower panel is at time t/ \JL/Ag ~ 4.2. Temperature contours are 
also shown on the faces of the simulation domain. At early times, 
conduction is more important than the RTI, and the tempera- 
ture contours on the faces of the simulation domain shown in the 
upper panel show that the temperature mixing layer has diffused 
substantially. The perturbation on the T = 38 surface are partly 
due to the tops of the RTI bubbles which are growing despite the 
diffusion. These bubbles soon emerge from the diffusion region, 
as seen in the lower panel. 
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Figure 12. The height of the temperature mixing layer (de- 
fined in *JsJ) as a function of time for our RTI-stable simula- 
tions with vertical initial magnetic fields. The solid lines show 
the results for simulations with anisotropic conductivity, whereas 
dashed lines show the results for simulations with isotropic con- 
ductivity. The top pair of lines, which are overlapping, are for 
Wc/wrti = 0.09 (VS2 and VS2I in Table [2j. The middle set of 
three lines all correspond to uj c /ui^i = 0.009. The solid line has 
anisotropic conductivity (VS3), the dashed line has isotropic con- 
ductivity (VS3I), and the dotted line has anisotropic conductivity, 
but is at higher resolution (VS3HR). The bottom two lines are 
for uv/wrti = 0.0009. For the simulations with relatively small 
anisotropic conductivities, the HBI becomes non-linear before the 
temperature mixing layer reaches the upper boundary of the sim- 
ulation. This results in deviations between the heights from the 
simulations with isotropic and anisotropic conductivity. The HBI 
is more effective at diminishing conduction in the high resolution 
simulation. 



5.2.1 RTI-stable, HBI-unstable (p+ < P-) 

We will first investigate the RTI-stable, HBI-unstable case. 
These runs simulate the discontinuous limit of the HBI. The 
presence of the HBI in our simulations can be inferred from 
the height of the temperature mixing layer as a function of 

The solid (dashed) lines corre- 



12 



time, as shown in Fig. 
spond to simulations with anisotropic (isotropic) conductiv- 
ity. The dotted line shows the results of our high resolution 
simulation (VS3HR in Table[2|. We find that hr increases as 
\/xt for the simulations with isotropic conductivity. For the 
relatively high conductivity of uj c /ujkti = 0.09, the height of 
the temperature mixing layer is the same for the simulation 
with anisotropic conductivity (VS2) as for the simulation 
with isotropic conductivity (VS2I). This is because the HBI 
has not had enough time to modify the conduction across 
the entire simulation domain before the temperature mix- 
ing layer hits the upper boundary. However, for simulations 
with cj c /ojrti = 0.009 or 0.0009 the HBI is able to partially 
inhibit conduction, slowing the growth of the height of the 
temperature mixing layer. 

The HBI saturates by reorienting vertical field lines to a 



more horizontal geometry (Parrish & Quataert 2008). Hor- 
izontal fields are produced by strong horizontal flows ( |Mc-| 
|Court et al.|201lj ) . This reduces the amount of vertical heat 
conduction, slowing the growth of the temperature mixing 
layer. In the rapid conduction limit (oj c <; i^buoy), the satu- 
ration time-scale is several times the buoyancy time-scale, 
^buoy = \fgd l°g T/dz. Conversely, if conduction is slow 
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magnetic field affects the evolution of hr- In the HBI non- 
linear evolution, the magnetic field orientation changes as 
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Figure 13. 01) as a function of height at various times for our 
high resolution simulation with cj c /(^>rti = 0.009 (VS3 in Table 
B. We plot (b 2 z ) at t/y/L/Ag = 28.2 (solid line), 56.6 (dashed 
line), and 84.8 (dot-dashed line). The horizontal dotted line shows 
(b z ) = f/3, which corresponds to an isotropic magnetic field. Due 
to the action of the HBI, the magnetic field becomes increasingly 
horizontal as the simulation progresses in time, especially near the 
centre of the domain where the buoyancy time is the shortest. 



(uj c < ujbuoy), the saturation time-scale is several times the 
conduction time-scale. The buoyancy frequency in our sim- 
ulations is initially very high because of the temperature 
discontinuity, so the HBI initially grows on the conduction 
time-scale. Thus, the HBI initially grows more slowly in runs 
with lower conductivity than in runs with higher conductiv- 
ity. This is why the HBI has had a more prominent effect 
on the height in the simulations with oj c /u;rti = 0.009 than 
in the simulations with o; c /o;rti = 0.0009 (see Fig. 121. We 



find that our high resolution simulation exhibits the slowest 
increase in height of all the HBI-unstable simulations. Due 
to increased resolution, the field is able to become more hor- 
izontal than in the simulations with our fiducial resolution, 
thereby decreasing the effective vertical conductivity. 

Previous simulations of the HBI show that although 
the field lines never become horizontal, the angle they form 



with the ho rizontal decreases as ~ (t/tbuoy) 
et al onl 1 " ~-t, + _ . . — i 



( McCourt 



20111, where ibuoy ~ ^buoy- To measure the average 
orientation of the magnetic field relative to the horizontal, 
we define (b 2 ) to be 

Bl 
BB 



(30) 



where (■) denotes horizontal average, as above. For simula- 
tions with initially vertical magnetic fields, (b 2 ) is initially 
equal to on e ev erywhere. 

In Fig. [13] we plot (bl) as a function of height for our 
high resolution simulation with o> c /u;rti = 0.009 (VS3HR 
m Table [2). If thermal conduction is isotropic, (b 2 ) remains 
close to one for the entire simulation. Figure [l3]shows, how- 
ever, that (b 2 ) decreases due to the HBI. We know that the 
decrease in (b 2 ) is not due to random mixing, as (b 2 ) de- 
creases to below the isotropic value of 1/3. There is also an 
asymmetry between (b 2 ) for z < and for z > 0. Because 



the pressure force pushes the interface upwards (see [ 5.1.1 1, 
the HBI acts somewhat more strongly for z > than for 
z < 0, causing (b 2 ) to be smaller for z > than for z < 0. 
We can estimate how the change in orientation of the 



(bl) 



\tj ~tbuoy) 



(31) 



As seen in Fig. 13 (6^) is a strong function of height. How- 
ever, for simplicity we will only consider (6J) at z = 
and only consider heat conduction through the z — 
plane. The energy equation implies dhr /dt ~ Xcff / hr , where 
Xoff = x("z}- We assume dlogT/dz ~ 1/hr and thus equate 
tbuoy ~ \fhr/g- One can then show 



dhr 
dt 



'X 



1 



(32) 



This implies that hr should converge to a finite value at 
long times, specifically 

y 2/3 

h T ~ fiTF- ( 33 ) 

In this calculation we have neglected terms of order unity. 
It is unclear if hr is actually converging to a finite value 
in our simulations that show clear evidence for the HBI. 
This is not surprising as even our simulations that start 
with h orizontal magnetic fields conduct heat vertically (see 
[ 5.1.1 1. Put another way, finite resolution prevents (b 2 ) from 
decreasing below some minimum value. 

In our RTI-stable but HBI-unstable simulations with 
initially vertical magnetic fields, we find that the evolution 
of the heights associated with our passive scalars, hc L and 
hc v , are qualitatively similar to the results from the RTI- 
stable simulations with horizontal initial magnetic fields. 
Furthermore, the heights are very similar when comparing 
simulations with isotropic and anisotropic conduction. This 
is expected at early times, as there is little mixing occurring 
at the interface until the HBI becomes non-linear. Even in 
the non-linear phase it is not that surprising that the HBI 
does not contribute significantly to mixing at the interface: 
as described in |McCourt et al.| ( |2011[ >, the HBI does not 
excite large velocities during its saturation. 

In Fig. [14] we plot isosurfaces of Cu for the high res- 
olution HBI-unstable simulation with w c /wrti = 0.009 
(VS3HR in Table (2). The vertical extent of the composi- 
tion mixing layer is much smaller than the height of the 
temperature mixing layer, which is about 0.6 (see Fig. |12[ ). 
There are also structures moving horizontally due to the flow 
produced by the HBI. This motion produces the horizontal 
magnetic fields which are the saturated state of the HBI. 

We can see the magnetic field growth associated with 
the HBI in Fig. [8] Our high resolution RTI-stable simula- 
tion with a vertical initial magnetic field (triple-dot-dashed 
line) experiences linear growth in magnetic energy starting 
at tj \/L/Ag » 7. By the end of the simulation the magnetic 
energy has increased by a factor of nine. The simulation with 
an initially vertical magnetic field but isotropic conduction 
exhibits essentially no magnetic field growth. Thus, this lin- 
ear magnetic field growth is due to the HBI. Because the 
HBI does not generate large velocities, the growth is not as 
extreme as in the RTI-unstable simulations. 



5.2.2 RTI-unstable, HBI-stable (p+ > p-J 

The RTI-unstable case with a vertical magnetic field is very 
similar to the RTI-unstable case with an initially horizontal 
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Figure 14. Isosurfaces of fixed composition Cjj at Cjj = 0.95 
and Cxj = 0.05 for the high resolution RTI-stable, HBI-unstable 
simulation with w c /w RTI = 0.009 (VS3HR in Tabled at time 
t/^/L/Ag = 71. Contours of Cjj are also shown on the faces of the 
simulation domain. For simulations with isotropic viscosity Cjj 
has very little horizontal variation. Thus, any horizontal structure 
is due to the HBI. 



magnetic field ({5.1.21. We plot the height of the temper- 
ature mixing layer as a function of time in Fig. |15[ The 
results are very similar to the corresponding plot in the 
RTI-unstable case with an initially horizontal magnetic field 
(Fig. [9). The results are similar because initially both sets of 
simulations are dominated by vertical heat conduction, but 
then transition to the RTI at later times. The most promi- 
nent difference is that in the initially vertical magnetic field 
case, when the RTI becomes important the height becomes 
very close to the height in the fiducial simulation. This is 
because vertical heat conduction does not impede the rise 
of RTI bubbles as much as horizontal heat conduction, so 
diffusive effects slow the rise of RTI bubbles less in the simu- 
lations with initially vertical magnetic fields. The RTI bub- 
bles in simulations with initially vertical fields are thinner 
than the RTI bubbles in simulations with initially horizontal 
fields, as in our simulation without explicit conductivity (see 
Fig.j5|. As in the RTI-unstable simulations with an initially 
horizontal magnetic field, the height in the high resolution 
simulation grows slightly more slowly than at the fiducial 
resolution. 

The magnetic field growth for the RTI-unstable case 
with initially vertical field is also similar to the case with ini- 
tially horizontal field (see Fig. [8]| . As for all the RTI-unstable 
simulations, the magnetic energy increases exponentially, 
and with about the same growth rate as for the simula- 
tions with initially horizontal magnetic fields. However, the 
exponential growth of magnetic energy begins much later 
(at tj \J Lj Ag « 3 rather than t ~ 0) in the simulation with 
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Figure 15. Heights and depths of the temperature mixing layer 
(defined in jj5| as a function of time for simulations with initially 
vertical magnetic fields, and with various cj c /ujrti- The solid lines 
show hx for our simulations which are RTI-unstable, whereas the 
dotted lines show — dj- for our simulations which are RTI-stable, 
and which have isotropic conductivity. The top pair of lines (in 
blue), which are overlapping, are for u> c /u)-g,Ti = 0.9 (VI and VS1I 
in Table[2j|. The next pair of lines (in red) are for Uc/k'RTI = 0.09 
(V2 and VS2I). The last set of two lines (in yellow) correspond to 
^c/^rti = 0.009 (V3 and VS3I). The bottom solid line shows hx 
for our fiducial simulation without explicit conductivity (RT in 
Table [TJ. The heights for the RTI-unstable runs are about equal 
to the depths for the RTI-stable runs when conduction dominates 
over the RTI. When the RTI becomes important, the heights in 
the RTI-unstable simulations increase in the same way as the 
height in the fiducial simulation. 



an initially vertical magnetic field, and the magnetic energy 
only amplifies by about 10 2 . This is because there is not as 
much horizontal motion to bend and amplify the initially 
vertical field lines. 



5.3 Skewed magnetic fields 

We have also run a set of simulations with initial magnetic 
fields which make a 45 degree angle with the horizontal (runs 
A3 and AS2 in Table[2} . These were compared to simulations 
with initially vertical magnetic fields, but with anisotropic 
thermal conduction with a magnitude one half as large (so 
the initial vertical heat conduction, proportional to x(^z}> 
was kept constant). The RTI-stable simulation with an ini- 
tially skewed magnetic field (AS2) was run with a conductiv- 
ity of tj c / 'ojrti = 0.18, and is very similar to the RTI-stable 
simulation with an initially vertical magnetic field, but with 
c^c/ujrti = 0.09 (VS2). We did not run the simulation with 
a sufficiently low conductivity to see the effects of the HBI, 
which are presumably not as strong for the simulation with 
initially skewed magnetic fields. 

For RTI-unstable simulations, the simulation with an 
initially skewed magnetic field and cj c /cjrti = 0.018 (A3), 
and the simulation with an initially vertical magnetic field 
and o; c /cjrti = 0.009 (V3) have the same height initially, 
when the temperature mixing layer is primarily expanding 
due to diffusion. However, when the RTI becomes dominant, 
the height for the simulation with an initially skewed mag- 
netic field is somewhat lower than the height for the simula- 
tion with the initially vertical magnetic field. This is because 
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in the RTI-dominated regime, the magnitude of the diffusiv- 
ity matters more than the direction of the initial field, as the 
magnetic field becomes isotropized. The simulation with an 
initially skewed magnetic field has uj c /lob,ti twice as large 
as the simulation with an initially vertical magnetic field, 
so this simulation is more dissipative. This causes a slower 
increase in height. 

In physical systems, the magnetic field will never be 
perfectly aligned or perpendicular to the contact disconti- 
nuity. However, these results show that if the magnetic field 
makes a moderate angle with the contact discontinuity, the 
behavior is very similar to a system with a vertical mag- 



netic field (see [ 5.2 I, but with a smaller conductivity. If the 



magnetic field makes a small angle with the contact discon- 
tinuity, then the evolution will likely be as described by the 
horizontal magnetic field simulations (see [5.1 1, which did 
have some vertical temperature diffusion. 



5.4 Additional physical effects 

Thermal conduction is anisotropic when the electron gyro- 
radius is much smaller than the electron mean free path. 
When this is true, the ion gyroradius is often also smaller 
than the ion mean free path. In this case, viscosity also acts 
predominantly along magnetic field lines. Because momen- 
tum is carried by the ions, the viscosity is smaller than the 
thermal diffusivity by a factor of about 50 or 100 if [i = 0.5 
or 1, respectively. We carried out simulations implementing 
anisotropic viscosity (as described in Parrish et al.|2012 1 for 
contact discontinuities where are RTI-unstablc. We found 
that the height of the temperature mixing layer is not sig- 
nificantly affected by inclusion of anisotropic viscosity. This 
is probably because the viscosity is so much smaller than the 
thermal diffusivity For simulations with an initially horizon- 
tal magnetic field, the RTI bubbles were more elongated in 
the direction of the initial magnetic field than in the sim- 
ulations with only anisotropic conductivity and numerical 
viscosity (see also |Dong fc Stone||2009| . 



6 SUMMARY AND DISCUSSION 

In a dilute magnetized plasma, thermal conduction is 
anisotropic; this can lead to qualitative changes in the be- 
havior of the plasma. We have investigated the evolution 
of RTI-stable and RTI-unstable contact discontinuities with 
anisotropic thermal conduction. The linear problem is only 
well-posed for a horizontal initial magnetic field. In this case, 
there is no change to the RTI dispersion relation in the local, 
Boussinesq limit. In this limit, the classical RTI perturba- 
tions are isothermal so they are unaffected by anisotropic 
conduction (eqn. [T8| ). More generally, compressible pertur- 
bations are affected by anisotropic conduction (eqn. 22 1. If 



the conduction time-scale is short compared to the growth 
time of the RTI, perturbations are isothermal and grow 
faster than the adiabatic perturbations of the classical RTI. 
The enhancement in the growth rate is modest, but is a 
function of the temperature contrast across the contact dis- 
continuity and the compressibility of the mode (see Fig. [2j|. 

We have run a number of fully non-linear simulations 
of contact discontinuities that are either RTI-stable or RTI- 
unstable, and with several initial magnetic field geometries 



(see Fig. [T]). We focus on initially horizontal and initially 
vertical magnetic fields - these two choices of the magnetic 
field orientation bracket the range of dynamics associated 
with contact discontinuities in dilute plasmas. We primarily 
used the heights of the temperature and composition mixing 
layers as diagnostics for the evolution of the initial temper- 
ature and density discontinuity. Interestingly, in our high 
resolution simulation with no explicit conduction (i.e., sim- 
ulating the classic RTI) the height increases linearly with 
time in the non-linear phase, rather than quadratically (see 
Fig. j3j|. This result is independent of the exact definition 
of height we use, or when we assume self-similar non-linear 
growth begins. 

In a dilute plasma, the non-linear evolution of a con- 
tact discontinuity can be described as a combination of three 
processes: vertical temperature diffusion, the RTI, and the 
HBI. For the simulations of RTI-stable contact discontinu- 
ities with anisotropic conduction, vertical temperature diffu- 
sion is most important. This is true even in our simulations 
with horizontal initial magnetic fields. The diffusion in this 
case occurs because we initially perturb the vertical velocity, 
and subsequent numerical mixing across the interface. Thus, 
both in simulations with vertical initial magnetic fields and 
horizontal initial magnetic fields, the height of the tempera- 
ture mixing layer increases due to diffusion. Although there 
is significant temperature diffusion, there is not much mix- 
ing at the interface in the absence of the RTI, so the height 
of the composition mixing layer remains small (see Fig. j7j|. 

An RTI-stable contact discontinuity is unstable to the 
HBI in the presence of vertical magnetic fields; the rate of 
vertical temperature diffusion across the contact disconti- 



nuity slows significantly due to the instability (see Fig. 12 1. 
The HBI reorients the vertical magnetic fields to a more 
horizontal geometry, impeding further heat conduction. We 
estimate that the HBI would cause a contact discontinuity 
to spread to a finite width ~ X 2 ^ /g 1 ^ ■ 

Our simulations with RTI-unstable contact discontinu- 
ities initially evolve in the same way as the simulations with 
RTI-stable contact discontinuities - at early times, the sim- 
ulations are dominated by temperature diffusion (see, e.g., 
Fig. IjjJ. However, RTI bubbles are forming and beginning 
to rise even as the temperature discontinuity diffuses away. 
The RTI bubbles rise faster and overtake the temperature 
diffusion, leading to dynamics dominated by the RTI at late 
times. At this stage, anisotropic conductivity does not play 
an especially important role. The heights of the RTI bubbles 
are similar in the simulations with and without anisotropic 
conductivity. However, the RTI bubbles are a little lower 
in the simulations with anisotropic conductivity, due to in- 
creased dissipation in the system. This is especially true for 
the simulations with initially horizontal magnetic fields. In 
this case, the RTI bubbles leak heat to their sides, becoming 
less buoyant. 

We now discuss astrophysical contexts in which our re- 
sults may be applicable, specifically, supernova remnants, 
emerging magnetic flux in the solar corona, and cold fronts 
in galaxy clusters. 

|Tilley fc Ba lsara ( 2006| discuss the effects of anisotropic 
conduction in a supernova remnant. They simulated an en- 
tire supernova remnant, but unlike Chevalier et al. ( 1992 1 



did not find any RTI. Surprisingly, it does not seem that 
there is significant temperature conduction across their con- 
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tact discontinuity, as the remnant maintains its spherical 
shape, instead of becoming more cylindrical as one would 
expect with their initial magnetic field. 

In supernova remnants, the contact discontinuity feels 
an effective gravitational force due to its deceleration into 
the ambient medium. The magnitude o f the accelera tion is 

~ 8/5 (e.g.,|Sedov|l946|. Thus, 



decreasing with time as g c R ~ t~ 
one might expect the height of RTI bubbles to increase as 
ffcfti 2 ~ t 2//5 . At late times this is slower than the diffusive 
growth of the height of the temperature mixing layer, which 
grows as */xt. As an example, assume a supernova remnant 
with an electron temperature of T — 10 8 K and number 
density n = 4 cm -3 is expanding into a medium with n 
I cm 3 . Then we have that the height of the temperature 
mixing layer due to temperature diffusion after a Sedov time 
~ 100 yr is 



\t ~ 3 pc 



t 



100 yr 



1/2 



(34) 



After ~ 100 yr the remnant would expand to a radius of 
about r ~ 2 pc, so g e g ~ r/t 2 ~ 0.6 cm s -2 . We find that 
the height of the highest RTI bubble after 100 yr would be 
about 



g e &t ~ 2 pc 



( t 



V 100 yr 



2/5 



(35) 



This suggests that diffusion and RTI dynamics are about 
equally important after a Sedov time, although diffusion 
would likely dominate at later times. Thus, a global sim- 
ulation of an entire supernova remnant would be necessary 
to accurately determine how anisotropic conduction affects 
the contact discontinuity. 



Berger et al. (20111 recently reported observations of 



hot, dilute solar prominences. They interpret the observa- 
tions as showing the RTI on the surface of the bubble of hot 
plasma, as described in numerical simulations (e.g., |Isobe 



|et al.|[2~005| ). They report temperatures of the hot plasma 
to be T fa 10 6 K, which, using n « 10 9 cm -3 , gives a col- 
lision frequency of v ~ 10 2 Hz. This is much smaller than 
the electron gyrofrequency uj cc « 10 7 (-B/1 G) Hz, indicating 
the heat conduction is primarily along magnetic field lines. 
Furthermore, such a plasma would have a parallel heat dif- 
fusivity of ~ 10 16 cm 2 s _1 (Spitzer 19621, so the diffusion 
time across a prominence which has a size » 10 9 cm is about 
100 s. This is much shorter than the evolution time for the 
prominence, which is on order an hour. Since conduction be- 
comes increasingly important at smaller length scales, this 
shows that anisotropic conduction should be important for 
the RTI discussed in |Berger et al.| \2Q11\ . 

The magnetic field in prominences is believed to be par- 
allel to the interface of the prominence. If some field lines 
were perpendicular to the interface, then heat conduction 
would quickly smear out the interface (see S 5.2.2 1. How- 
ever, the observations show that the prominences hold their 
shape on time-scales of an hour. Thus, we can take the mag- 
netic fields to be primarily along the interface. Our linear 
results in ij2] show that the growth of the RTI should be 
faster than would be predicted without anisotropic conduc- 
tion, though how much faster depends on the compressibility 



of the modes. In the non-linear regime ({ 5.1.2 1, we expect 
a superposition of spreading of the interface due to perpen- 
dicular heat diffusion and a somewhat slower RTI growth. 



Because perpendicular diffusion is not very important (the 
perpendicular diffusivity is smaller by a factor of (v/ui cc ) 2 ), 
we expect the non-linear behavior to be similar to the RTI 
without conduction. 

Our calculations are not ideally suited for compari- 
son with solar prominences as we have assumed that the 
magnetic fields are dynamically unimportant, which is not 
the case for solar prominences. Strong magnetic fields in- 
hibit motions which bend the field lines, instead favoring 
interchange-like motions (SG07). Anisotropic viscosity along 



field lines also would inhibit motion along field lines (Dong 
fc Stone|2009 i 



Parrish fc Quataert (2008) and Birnboim et al. (20101 



suggested that thermal diffusion across the contact discon- 
tinuities found at cold fronts in galaxy clusters could be 
suppressed by the HBI. In §5.2. 1| we determined that the 
HBI will cause the conta ct discontinuity t o diffu se to a fi- 



y 2/3 ,1/3 



Birnboim et al. 



(2010) simulate 



nite width 

shocks propagating through a galaxy cluster which merge 
producing a cold front. Using their cluster parameters for a 
front at 200 kpc, we find that x ~ 3 x 10 30 cm 2 s~\ and 
g w 10 -8 cm s -2 . This gives the saturation width of the 
contact discontinuity to be ~ 30 kpc, which is somewhat 
larger than the Chandra upper limit of 10 kpc. The satura- 
tion width of the cold front observed in the galaxy cluster 



A3667 (Vikhlinin et al. 20011 can be estimated to also be 



~ 40 kpc - much larger than the observed width ^ 5 kpc. 
Furthermore, our numerical results in §5.2.1| suggests that 
the saturation width is likely a few times \ ' 3 /g 1 ^ 3 . These 
results suggest that additional effects not included in our 
simple model might be required to explain the narrow widths 
of observed cold fronts. 



APPENDIX A: CALIBRATION OF LINEAR 
GROWTH RATES USING NUMERICAL 
SIMULATIONS 

We have tried to confirm the growth rates derived in Sj2] 
through numerical simulations. Our numerical methodology 
is summarized in f^3j There are two difficulties that make 
accurate comparison between linear theory and the numeri- 
cal results impossible. The first difficulty is that the fastest 
growing modes are those on the smallest length scale allowed 
in the simulation. The second difficulty is that viscous effects 
will act to slow the growth of the instability. Although we 
do not implement an explicit viscosity in our numerical sim- 
ulations (see ^3|, there is numerical diffusion which adds an 
effective viscosity to the problem. This viscosity will damp 
out any modes with wavelengths smaller than a cutoff wave- 
length, which is the width of a few grid points - we will re- 
fer to this cutoff wavelength as the viscous length scale. The 
fastest growing modes will be those with a wavelength only 
slightly larger than this viscous length scale. Thus, these 
modes are substantially affected by viscosity, and cannot be 
accurately compared to our calculations above which are for 
the inviscid RTI problem. 

A possible way to get around these difficulties would 
be to only excite long wavelength modes which are not af- 
fected by viscosity. However, rounding errors produce ran- 
dom perturbations at the smallest scales in the simulation, 
and these random perturbations grow faster than the long 
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wavelength modes of interest (as their growth rates are many 
times larger than the growth rates of the modes of interest). 
If we increase the amplitude of the initial perturbations, it 
will take longer for the random noise to grow to the same size 
as the initial perturbation. Unfortunately, non-linear effects 
become important at these higher amplitudes, so we still 
cannot accurately measure the linear growth rate. 

We have tried using the eigenfunctions from the calcula- 
tion in fj2] as an initial perturbation, in hopes of seeing some 
growth due to the initial perturbation before the simulation 
is overwhelmed by the small scale modes. Unfortunately, the 
cusp in the perturbation vertical velocity and discontinuities 
in perturbation horizontal velocity, density, and pressure at 
z = are not resolved by the simulation. Because of this, 
using the eigenfunctions as an initial perturbation launches 
sound waves with amplitudes similar to the amplitude of the 
initial perturbation. These obfuscate the growth of the mode 
of interest. To summarize, we have sound waves launched by 
the initial perturbation which dominate the dynamics at the 
beginning of the simulation, and quickly growing small scale 
modes which dominate later on in the simulation. Between 
these two effects, we are unable to measure the growth of the 
longer wavelength modes which are relatively unaffected by 
viscosity. Presumably these difficulties could be remedied 
by using a resolved density profile rather than a true dis- 
continuity on the grid scale. However, we concluded that 
there was not sufficient motivation to do so given that our 
numerical methods have been validated using many other 



tests (in particular, see Parrish & Stone ( 2005 I and Sharma 
& Hammett (20071 for tests of the anisotropic conduction 



methods). 
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